#include "func_EHw.h"
#include <gsl/gsl_math.h>
#include <gsl/gsl_const_mksa.h>




struct subint_context {
	enum field_EH field_class; // E or H
	uint8_t j; // component: 0 - x, 1 - y, 2 - z
	double w; // frequency
	double r[3]; // Observation point
	
	const struct current * crnt;
	const union trajectory * trj;
	const union trajectory_len * trlen;
};
/**
 * Vector conjugation. 
 * \ksi=[\zeta,\eta]
 * \ksi_i = \epsilon_{ikl}\zeta_k\eta_l
 * @param i - component to output
 * @param zeta - first argument (double[3])
 * @param eta - second argument (double[3])
 * @return vector conjugate
 */
static double vconj(uint8_t i, const double * zeta, const double * eta){
	switch(i){
		case 0:
			return zeta[1]*eta[2] - zeta[2]*eta[1];
		case 1:
			return -zeta[0]*eta[2] + zeta[2]*eta[0];
		case 2:
			return -zeta[1]*eta[0] + zeta[0]*eta[1];
		default:
			return GSL_NAN;
	}
}

/**
 * Calculates the expression to be integrated over. Fucntion arguments
 * match the specification from func_trapz.h
 * @param t - parameter \in [0,1]
 * @param context - pointer to subint_context
 * @return complex value of the expression, being calculated
 */
static complex subint_func(double t, void * ctx){
	struct subint_context * data = ctx;
	double w = data -> w;
	double k = w / GSL_CONST_MKSA_SPEED_OF_LIGHT; // We are in ~ air
	double epsilon = GSL_CONST_MKSA_VACUUM_PERMITTIVITY; // We are again in ~ air
	double l = trajectory_len_calc(data -> trlen, t);
	complex I_w = current_get_spectrum(data -> crnt,w,l);
	double R_tau[3]; // Current trajectory point coordinates
	trajectory_calc(data -> trj, t, R_tau, TRAJECTORY_CALC_R);
	double Rv[3];//vector R
	uint8_t j;//A nice index to use in for cycles
	for(j=0;j<3;j+=1){
		Rv[j] = data -> r[j] - R_tau[j];
	}
	double n[3];
	double R = gsl_hypot3(Rv[0], Rv[1], Rv[2]);
	for(j=0;j<3;j+=1){
		n[j] = Rv[j]/R;
	}
	double dr_dt[3], dl_dt; // needed for using t as integration variable instead of l
	trajectory_calc(data -> trj, t, dr_dt, TRAJECTORY_CALC_DR_DT);
	dl_dt = gsl_hypot3(dr_dt[0], dr_dt[1], dr_dt[2]);
	
	double tau[3];
	for(j = 0; j < 3; j+=1){
		tau[j] = dr_dt[j] / dl_dt; // Normalizing derrivative vector to tau
	}
	
	j = data -> j; // Now j is the projection to be got
	complex ans;
	switch(data -> field_class){
		case FIELD_E:{
			complex E_jw = I*k/(4*M_PI*epsilon*w)*
				I_w * (/*[*/  (k * cexp(I*k*R))/R*(tau[j] - n[j]*sc(tau,n)) + 
				I*cexp(I*k*R)/pow(R,2)*(tau[j] - 3*n[j]*sc(tau,n))*(1+I/(k*R)) /*]*/);
			ans = dl_dt * E_jw;
		}
		break;
		case FIELD_H:{
			complex H_jw = 1/(4*M_PI) * I_w * vconj(j,n,tau) * cexp(I*k*R)/(R)*(I*k - 1/R);
			ans = dl_dt * H_jw;
		}
		break;
		default:
			ans = GSL_NAN;
	}
	return ans;
}



int func_EHw_calc(enum field_EH field_class, double eps, const double * r, const double * n,
	double w, const union trajectory * trj, const union trajectory_len * trlen, 
	const struct current * crnt, complex * result){
	
	// Which components will we need? Yet it's OK to check n_i == 0.000, 
	// because no floating point errors have yet accumulated, that would
	// spoil the comparisson.
	
	uint8_t calculate_projections = 0;
	uint8_t j;
	for(j=0;j<3;j+=1){
		if(n[j] == 0.000){
			//Bit allready zero
		}else{
			calculate_projections|= 1 << j; //Set bit j to '1'
		}
	}
	
	// Now let's normalize n
	double n_n = gsl_hypot3(n[0], n[1], n[2]);
	double nn[3]; // n_n - n norm, nn - normalized n
	for(j=0;j<3;j+=1){
		nn[j] = n[j] / n_n;
	}
	// From now we should use nn, not n
	
	struct subint_context ctx;
	ctx.field_class = field_class;
	ctx.w = w;
	copyvec(ctx.r, r);
	ctx.trj = trj;
	ctx.trlen = trlen;
	ctx.crnt = crnt;
	
	complex sum = 0; // For all active projections the sum
	for(j=0;j<3;j+=1){//For all projections
		if( ! ((calculate_projections) & (1 << j)) ){
			continue;// Skip projection
		}
		ctx.j = j;
		complex intval;
		int integration_status = 
			func_trapz(subint_func, &ctx, 0, 1, eps, &intval, NULL, NULL);
		if(integration_status){
			* result = GSL_NAN;
			return integration_status;
		}
		// Integrations seems OK
		sum += nn[j] * intval; // Scalar (n,field{E,H}_w)
	}
	// OK!
	* result = sum;
	return 0;
}
